# Load packages
library(tidyverse)
library(sf)
# devtools::install_github("mabecker89/leaflet@crosstalk6")
library(leaflet)
library(leaflet.extras)
library(crosstalk)
library(DT)
library(sp)

# Set options for map size 
knitr::opts_chunk$set(fig.width = 12, fig.height = 10)

# Import data
cwpt_qs_all <- st_read("./Data/Spatial/Clean/cwpt_qs_all.shp",
                       stringsAsFactors = FALSE, quiet = TRUE)

cwpt_qs_mun_w50 <- read_csv("./Data/Clean/cwpt_qs_mun_w50.csv")

ab_counties <- st_read("./Data/Spatial/Clean/ab_counties_md_sa.shp",
                       stringsAsFactors = FALSE, quiet = TRUE)

# Set AOI
aoi_name <- "Lac Ste. Anne County"

# Filter for AOI in QS data
qs_aoi <- cwpt_qs_mun_w50 %>%
  filter(County_MD_SA == aoi_name) %>%
  select(County_MD_SA, LLD, 
         Quarter, Section, Township, Range, Meridian,
         j_land_mean, j_land_max, q_land_mean, q_land_max,
         j_pix_p_mean, q_pix_p_mean, 
         pct_rp:pct_ma, pct_ng,
         LinkID) %>%
  mutate(total_ag = pct_rp + pct_crop + pct_tp + pct_ma) %>%
  mutate(total_ag = if_else(total_ag > 1, 1, total_ag))

# Join QS data to QS geometry
cwpt_qs_aoi <- cwpt_qs_all %>%
  filter(LinkID %in% qs_aoi$LinkID) %>%
  left_join(qs_aoi, by = "LinkID") %>%
  st_transform("+init=epsg:4326")

# Filter for AOI in county shapefile
aoi <- ab_counties %>%
  filter(C_MD_SA == aoi_name) %>%
  st_transform("+init=epsg:4326")

names(st_geometry(cwpt_qs_aoi)) = NULL
names(st_geometry(aoi)) = NULL
# Colours
land <- colorNumeric("Blues", domain = cwpt_qs_aoi$j_land_mean)

# Generate leaflet map of landscape importance score
map_land_aoi <-
  aoi %>%
  leaflet() %>%
  addTiles() %>%
  addProviderTiles("Esri.WorldImagery", group = "Imagery") %>%
  addFullscreenControl() %>%
  addResetMapButton() %>%
  addMapPane(name = "AOI Boundary", zIndex = 420) %>%
  addMapPane(name = "Layers", zIndex = 410) %>%
  addScaleBar(position = "bottomleft",
              options = scaleBarOptions(imperial = FALSE)) %>%
  
  # Polygon layers
  addPolylines(color = "#070707", weight = 3, smoothFactor = 0.2,
               options = leafletOptions(pane = "AOI Boundary")) %>%
  
  addPolygons(data = cwpt_qs_aoi, color = "#444444", weight = 0.4, 
              smoothFactor = 0.5, opacity = 1.0, fillOpacity = 0.6,
              # Using UPSLOPE MEAN & JENKS
              fillColor = ~ land(j_land_mean),
              group = "Landscape Position Importance Score",
              popup = paste("Legal Land Description:", 
                            "<b>", cwpt_qs_aoi$LLD, "</b>", "<br>",
                            "Landscape Position Importance Score:",
                            "<b>", cwpt_qs_aoi$j_land_mean, "</b>"),
              highlightOptions = highlightOptions(color = "white", weight = 2,
                                                  bringToFront = TRUE),
              options = leafletOptions(pane = "Layers")) %>%
  
  # Legend layer
  addLegend(data = cwpt_qs_aoi, position = "bottomright", pal = land, 
            values = ~ j_land_mean, title = "Index Score", opacity = 1,
            group = "Landscape Position Importance Score") %>%
  
  # Layers control
  addLayersControl(overlayGroups = c("Landscape Position Importance Score",
                                     "Imagery"),
                   options = layersControlOptions(collapsed = FALSE)) %>%
  
  hideGroup("Imagery")

map_land_aoi
# Colours
cp <- colorNumeric(c("#00441b","#238b45", "#74c476","#a1d99b", "#e5f5e0", 
                     "#f6e8c3","#dfc27d", "#bf812d","#8c510a", "#543005"),
                   domain = cwpt_qs_aoi$j_pix_p_mean)

# Convert sf object to sp (required for crosstalk reactivity)
cwpt_qs_aoi_sp <- cwpt_qs_aoi %>%
  as_Spatial()

# Create SharedData objects
sd <- SharedData$new(cwpt_qs_aoi_sp)
sd_df <- SharedData$new(as.data.frame(cwpt_qs_aoi_sp@data), group = sd$groupName())

# Create filter slider
filter_slider("pixelscore", "Filter for Landscape Position Importance Score", 
              sd_df, ~j_land_mean, step = 0.5, width = 550, ticks = TRUE)
filter_slider("ag", 
              "Filter for proportion of quarter-section classified as Agriculture",
              sd_df, ~total_ag, step = 0.05, width = 550, ticks = TRUE)
# Generate leaflet map of on-pixel scores
map_pix_aoi <-
  leaflet(sd) %>%
  addTiles() %>%
  addProviderTiles("Esri.WorldImagery", group = "Imagery") %>%
  addFullscreenControl() %>%
  addResetMapButton() %>%
  addScaleBar(position = "bottomleft",
              options = scaleBarOptions(imperial = FALSE)) %>%
  addMapPane(name = "AOI Boundary", zIndex = 420) %>%
  addMapPane(name = "Layers", zIndex = 410) %>%
  
  # Add polygons
  addPolygons(color = "#444444", weight = 0.4, smoothFactor = 0.5,
              opacity = 1.0, fillOpacity = 0.6,
              # Using JENKS 
              fillColor = ~cp(j_pix_p_mean),
              group = "On-Pixel Contribution",
              popup = paste("Legal Land Description:", 
                            "<b>", cwpt_qs_aoi_sp$LLD, "</b>", "<br>",
                            "Pixel Contribution Score:",
                            "<b>", cwpt_qs_aoi_sp$j_pix_p_mean, "</b>"),
              highlightOptions = highlightOptions(color = "white", weight = 2,
                                                  bringToFront = TRUE),
              options = leafletOptions(pane = "Layers")) %>%
  
  addPolylines(data = aoi, color = "#070707", weight = 2.5, smoothFactor = 0.2,
               options = leafletOptions(pane = "AOI Boundary")) %>%
  
  # Add legend
  addLegend(position = "bottomright", pal = cp, values = ~ j_pix_p_mean,
            title = "Pixel Contribution Score", opacity = 1) %>%
  
  # Layers control
  addLayersControl(overlayGroups = c("Imagery", "On-Pixel Contribution"),
                   options = layersControlOptions(collapsed = FALSE)) %>%
  hideGroup("Imagery")
  
map_pix_aoi
datatable(data = sd_df,
          extensions = "Buttons",
          filter = "top",
          rownames = FALSE,
          caption = NULL,
          class = 'compact',
          width = "100%",
          options = list(pageLength = 10,
                         columnDefs = list(list(visible = FALSE,
                                              targets = c(0,2,9:11,13,17,19))),
                         dom = 'lBfrtip',
                         buttons = 
                           list(list(extend = 'collection',
                                     buttons =
                            list(list(extend = 'csv',
                                      filename = 'cwpt_abmi'),
                                 list(extend = 'excel',
                                      filename = 'cwpt_abmi')),
                         text = 'Download data'))),
          colnames = c("", "Municipality", "", 
                       "Quarter", "Section", "Township", "Range", "Meridian",
                       "Landscape Importance Score", "", "", "",
                       "On-Pixel Score", "",
                       "Rough Pasture",
                       "Cropland",
                       "Tame Pasture",
                       "",
                       "Native Grassland", "")
          ) %>%
  formatPercentage(c("pct_rp",
                     "pct_crop",
                     "pct_tp",
                     "pct_ma",
                     "pct_ng",
                     "total_ag")) %>%
  formatStyle(columns = c(4,8), textAlign = 'right')